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We study the phase diagram of the one-dimensional Hubbard model with next-nearest-neighbor 
hopping using exact diagonalization, the density-matrix renormalization group, the Edwards vari- 
ational ansatz, and an adaptation of weak-coupling calculations on the two-chain Hubbard model. 
We find that a substantial region of the strong-coupling phase diagram is ferromagnetic, and that 
three physically different limiting cases are connected in one ferromagnetic phase. At a point in 
the phase diagram at which there are two Fermi points at weak coupling, we study carefully the 
phase transition from the paramagnetic state to the fully polarized one as a function of the on-site 
Coulomb repulsion. We present evidence that the transition is second order and determine the crit- 
ical exponents numerically. In this parameter regime, the system can be described as a Luttinger 
liquid at weak coupling. We extract the Luttinger-liquid parameters and show how their behavior 
differs from that of the nearest-neighbor Hubbard model. The general weak-coupling phase dia- 
gram can be mapped onto that of the two-chain Hubbard model. We exhibit explicitly the adapted 
phase diagram and determine its validity by numerically calculating spin and charge gaps using the 
density-matrix renormalization group. 



PACS Numbers: 71.10.Fd, 75.10.Lp and 75.40. Mg 
I. INTRODUCTION 

Unlike strongly correlated phenomena such as antifer- 
romagnetism or superconductivity which can be treated 
starting from a weak-coupling point of view, metallic fer- 
romagnetism is an intrinsically intermediate or strong 
coupling phenomenon. Because of this, the origin of 
metallic ferromagnetism is still poorly understood, even 
after decades of research. The simplest model of corre- 
lated electrons, the Hubbard model, was introduced sk 
multaneously by Gutzwillerjd Hubbarda and Kanamoritl 
in 1963 in order to study ferromagnetism. Indeed, at a 
mean-field level, the Hubbard model seems to be a good 
starting point for ferromagnetism, because the Stoner cri- 
terion predicts a ferromagnetic ground state for a wide 
range of parameters. However, the inclusion of correla- 
tion effects makes the conditions for the appearance of 
ferromagnetism much more restrictiveJ3 

There is one limit of the Hubbard model in which a fer- 
romagnetic state can be obtained within an exact treat- 
ment. For U = oo the ground state of the half-filled 
system has macroscopic degeneracy since all states with 
different spin S have the same energy. For bipartite lat- 
tices (such as the hypercubic and bcc lattices) in dimen- 
sion d > 2, and for fee and hep lattices with negative 
hopping integrals, Nagaokaa proved that when one hole 
is then added to the i-site system, the ground state of 
the model has maximum spin. It has not yet been proven 
possible to generalize Nagaoka's proof for the stability of 
ferromagnetism to a finite density of holes S ^ 1 — n > 0, 
making the treatment of the thermodynamic limit prob- 
lematic. One can, however, use the opposite approach 
and try to show that the ferromagnetic state is unstable 
by applying a suitable variational wave function. Start- 
ing from the fully polarized state, one can flip a spin in 



an appropriate manner and then see if the corresponding 
energy is lower than the ferromagnetic one. If a lower en- 
ergy can be found, the fully polarized ferromagnetic state 
is proven to be unstable. Recently, by choosing a highly; 
sophisticated variational state, Wurth and co-workersO 
were able to bring the critical hole density Sc above 
which the Nagaoka state is unstable at C/ = cx) down, 
to 6c = 0.25J| for a square lattice. Exact diagonalizatioro 
and DMRGQ studies suggest that Sc could be even lower. 
Hence it cannot be ruled out that Sc — for bipartite 
lattices, as is the. case for the hyper-cubic lattice in infi- 
nite dimensions.u In contrast, for non-bipartite lattices, a 
partly polarized mound state (ferrimagnetism) has been 
obtained by LiebBI 

For one-dimensional systems, the situation is ep^n less 
favorable for ferromagnetism. Lieb and MattialJ have 
proven that the ground state is unmagnetized for any real 
and particle-symmetric but otherwise arbitrary interac- 
tion. This theorem applies to a single band in d — 1, 
provided that the hopping is only between nearest neigh- 
bors and the interaction involves only densities. Since 
both conditions are fulfilled in the Hubbard model, its 
ground state in d = 1 cannot be ferromagnetic. 

In principle, the Hubbard model is obtained from an 
extreme truncation of a more general Hamiltonian de- 
scribing interacting electrons in a solid. Only the on- 
site interaction and one relevant band are kept. Nearest- 
neighbor interaction (e.g. direct exchange), band degen- 
eracy and the associated Hund's rule couplings are totally 
neglected. In addition, the non-interacting band struc- 
ture and density of states can be strongly affected by 
the lattice structure. In order to enhance the possibility 
of ferromagnetism, one can modify the simple Hubbard 
model by putting some of these neglected features back 
in. 



Strack and VoUhardttll have studied a Hubbard model 
to which they have added all possible nearest-neighbor 
interactions: the usual interaction between charge {V), 
the density-dependent hopping (X), the Heisenberg ex- 
change (F) and pair hopping {F'). They show that for a 
particular range of parameter values the model has max- 
imum total spin at half-filling. These arguments can also 
be extended to the Nagaoka case of one hole.E3 

Another option is to take a multi-band Hamiltonian 
with d-band degeneracy together with a Hund'ajule cou- 
pling between the different d orbitals. Okabell3 has in- 
vestigated the stability of the ferromagnetic state of such 
a model variationally, while Fleck and coworkers have 
studied a-similar model including next-nearest-neighbor 
hopping.Eil These authors claim that the Hund's rule cou- 
pling is necessary to obtain ferramagnetism. Very re- 
cently, Biinemman and coworkeralj have studied a two- 
band Hubbard model with a multi-band Gutzwiller wave 
function. They found that a ferromagnetic transition oc- 
curs at large interaction and stress that a finite value of 
the exchange interaction is also required. We will see here 
that ferromagnetism can be obtained in a non-orbitally- 
degenerate model. 

Mielkeu has proven the following theorem, equivalent 
to Hund's rule, for a Hubbard model with a flat band. 
If the model has an Af-fold degenerate single-particle 
ground state, then for any number of electrons N < M 
the fully polarized state (with total spin S* = -y) is a 
ground state of the system. Additional conditions that 
determine whether this ground state is unique are also 
given, f^liis theorem has also been extended to nearly flat 
bands.E3 Recently, Tasakitj has considered a two-band 
Hubbard model with next-nearest-neighbor hopping. He 
has proven that the quarter-filled system (average elec- 
tron density n = 0.5) is ferromagnetic for la^se enough 
on-site interaction U. Penc and coworkersta have ex- 
tended this result to other fillings by studying a Hamil- 
tonian in which a chemical potential is added to all even 
sites of the lattice to make a perturbative argument valid. 
We shall see below that this term is unnecessary to obtain 
a ferromagnetic ground state. 

Several authors have studied the ordiiiary Hubbard 
model on various lattice structures. UlmkeEj investigated 
the case of an fee lattice using dynamical mean-field the- 
ory and quantum Monte Carlo simulations. He found 
ferromagnetism for intermediate values of U, using the 
density of states of the infinite-dimensional system. For 
the three-dimensional density of states he found that one 
must add next-nearest-neighbor hopping in order to ob- 
tain ferromagnetism. His general conclusion is that a nec- 
essary condition for ferromagnetism is a density of states 
with large spectral weight near the lower band edge. 
Hanisch and coworkersEil have investigated the stability 
of saturated ferromagnetism using a variational approach 
for various lattice structures in two and three dimen- 
sions. Their conclusions are similar to the previous ones, 
namely that a particle-hole asymmetry and a divergent 
density of states at the lower band energy are necessary 



ingredients for obtaining a ferromagnetic ground state. 
Similar conclusions have been reached very recently us- 
ing the spectral density approachE3 and dynamical iwan- 
field theory with the Non-Crossing Approximation.L3 

In the small density limit, the ferromagnetic state of 
the Hubbard model with arbitrary non-diagonal hop- 
ping and with a band structure with a quadratic dis- 
persion about the band minima has been shown to be 
unstable for d > 3, while for d = 2 a small window 
of parameters for which tiie fully polarized state is not 
destabilized still remains .cJ Ipdeed, a projector quan- 
tum Monte Carlo calculationES yields ferromagnetism 
precisely in this allowed region for the two-dimensional 
Hubbard model on a square lattice with next-nearest- 
neighbor hopping. Fin|al|ly, a renormalization group cal- 
culation for this modeled also yields ferromagnetism in a 
particular regimer— ■ 

In this work,EZl we study perhaps the simplest 
case of a Hubbard model exhibiting ferromagnetism: 
the one-dimensional Hubbard model with an addi- 
tional next-nearest-|iieighbor hopiping. Previously, ex- 
act diagonalization,Ej variationalEa and Densit¥-Matrix 
Renormalization Group (DMRG) calculation^ on this 
model have already concluded that there is an exten- 
sive ferromagnetic phase, for large enough coupling U. A 
weak-coupling analysisEHI applied to this one-dimensional 
model leads to a phase with a spin gap which is the 
one-dimensional analog of a superconductor. Projector 
QMC and—DMRG calculations for the special case of 
half-fiUingEl have recently been carried out at weak to 
intermediate coupling and are consistent with the weak- 
coupling analysis. Here we treat the strong and weak- 
coupling phase diagrams of the model comprehensively 
with numerical and variational techniques, and link them 
by studying the phase transition using exact diagonaliza- 
tion and DMRG calculations. 

This paper is organized as follows. In Sec. ||, we dis- 
cuss the basic properties of the model, motivate the ex- 
istence of ferromagnetism by discussing calculations on 
three and four site clusters, and discuss exactly treatable 
limiting cases at interaction U — oo, some of which yield 
ferromagnetic ground states. The results of a variational 
calculation of the U = oo phase diagram using the Ed- 
wards variational ansatz are presented in Sec. III. While 



this technique was previously applied to this model in 
Ref. ^ at finite C/, here we discuss the U = oo phase 
diagram and include a determination of the total spin 
of the variational state. Exact diagonalization calcula- 
tions, presented in Sec. tt^, are used to illustrate the 
determination of the critical interaction strength Uc at 
the ferromagnetic transition, to determine the order of 
the phase transition, and to examine the scaling of Uc 
for small next-nearest-neighbor hopping t2- In Sec. M we 
use the DMRG to determine Uc as a function of den- 
sity for three different t2 regimes, discuss the behavior of 
the Luttinger liquid parameters at the transition from a 
Luttinger liquid regime to the ferromagnetic regime, and 
discuss the behavior of the spin-spin correlation function 



near the transition. We present a determination of the 
strong-couphng (U = cx)) and weak-couphng phase di- 
agrams in Sec. VI. The strong-couphng phase diagram 
is calculated using the DMRG, and the weak-coupling 
phase diagram is determined by adapting the results of 
Balents and Fishertj for the two-chain|JjIubbard model 
to this model, as suggested by FabrizioE3 A DMRG cal- 
culation of the spin and charge gaps at weak but finite 
U is then used to check the validity of the weak-coupling 
phase diagram. The DMRG calculations of Uc and the 
U = oo phase diagram were reported in a preliminary 
form in Ref. U% but here they are carried out with more 
accuracy and the discussion is extended. 



II. THE ii - fa CHAIN 

A. The model 

We consider the one-dimensional Hubbard model with 
next-nearest-neighbor hopping (see Fig. ||) with the 
Hamiltonian 



H ^-tiY^ (cl+i^Cj^ -I- h.c. 



For U — and periodic boundary conditions H can be 
diagonalized by Fourier transformation yielding 



H = y2e(k)cl 



fc,(T 



Cko 



with k an integer multiple of ^ and 

e{k) — —2ti cos k — 2^2 cos 2k. 



(2) 



(3) 



(1) 



The dispersion e(fc) will have one minimum at fc = for 
t2 > —0.25 and two minima at finite k for t2 < —0.25. 
Similarly, there will be maxima in e(fc) at fc = ±7r for 
^2 < 0.25, but will shifted away from k — tt for t2 > 0.25. 
Therefore, for small 1^2 1, e(fc) does not differ qualitatively 
from the ^2 = band structure, and there are two Fermi 
points for arbitrary electron density n. On the other 
hand, for ^2 < —0.25 and sufficiently small densities or 
for ^2 > 0.25 and sufficiently large densities, the Fermi 
surface has four Fermi points, namely ifc^i and ifc^a 
and, as we will see, can be mapped to a two-band model 
at weak coupling. The resulting U — ground-state 
phase diagram is depicted in Fig. ||. Note that we take 
the horizontal axis to be — ^2 in this and all subsequent 
phase diagrams in order to better display the ^2 < 0, 
n < 1 region on which we will primarily concentrate in 
this work. 




We will call this model the ii — f2 chain. 

-t, 
-t, 



FIG. 1. The fi - f2 chain. 

The summation is over all L sites and spin a. Here we 
will always take U positive and set h = 1. The sign of 
ti is arbitrary since a gauge transformation Cj — > e^^'-^Cj 
can reverse it, so we set ii = 1 without loss of generality, 
and measure all energies in units of ii. This Hamilto- 
nian conserves the number of particles, the total spin 
S and its projection onto the quantization axis, Sz- If 
a particle-hole transformation is applied to the system, 
the transformation ^2 ^ ~^2 is necessary to recover the 
original Hamiltonian. Therefore, the parameter regime 
n > 1, ^2 > maps to n < 1, ^2 < and the regime 
n > 1, ^2 < maps to n < 1, ^2 > (see Fig. g). Because 
a definite order of the particles is oq longer enforced when 
t2 ^ 0, the Lieb-Mattis theoren£3 does not apply and, 
indeed, we will see that we do find ferromagnetism. 

The effect of the sign of f2 can be understood by con- 
sidering the Hamiltonian of Eq. (1^) with L = ?>^ N — 2 
particles, and open boundary conditions. This three-site 
model, treated by Tasaki in Ref. p3^, has a ferromagnetic 
ground state only when t2 < 0. For this choice of sign, 
the triangular structure frustrates the antiferromagnetic 
order normally found in the Hubbard model. 




FIG. 2. U — phase diagram. The inset plots show the 
qualitative behavior of e{k) in the different regions with the 
Fermi level indicated by the dashed line. 



B. Square cluster 

It is useful to solve the Hamiltonian of Eq. (|l|) ex- 
actly for L = A and periodic boundary conditions in anal- 
ogy with the solution of the three-site problem discussed 
above. This will allow us to further examine the effect of 
the sign of ^2 on the ground state, and also will illustrate 
finite-size effects due to open and closed shell configura- 
tions. As depicted in Fig. H, one obtains the Hubbard 
model on a square with additional diagonal hopping. 




FIG. 3. Square cluster. 

We consider the U = oo limit with N = 2 and A^ = 3 
electrons and minimum Sz (e.g. Sz — for A^ = 2). 
The dimension of the Hilbert space is 12 in both cases. 
With the aid of group theory, the problem can be solved 
analytically. The Hamiltonian has the symmetry of the 
group Civ (see Table g), and can be diagonahzed using 
the symmetry-adapted wave functions. The eigenvalues 
are labeled by the corresponding irreducible representa- 
tions {Ai ,A2,Bi,B2,E). In order to see whether a state 
is fully polarized or not, we check whether the eigen- 
function is symmetric or antisymmetric. The irreducible 
representations.^ 1, _Bi and B2 are symmetric and A2 is 
antisymmetries Since the global wave function is re- 
quired to be antisymmetric and the spin function of a 
fully polarized state is symmetric, the ferromagnetic state 
must belong to the representation A2. 

An analysis of the eigenvalues (see Tables || and [D ) 
leads to the following conclusions: 

A?^ = 2 : For t2 > -0.5 the ground state has 5* = 0. In 
this case the fully polarized state is not a closed shell 
state and it seems that for this reason it can not be the 
ground state. Similar open shell effects have also been 
observed in numerical exact diagonalipation calculations 
for larger systems and other fillings.Ea 

For t2 < —0.5 (where the single-particle spectrum has 
two minima) the ferromagnetic {S = 1) ground state of 
representation A2 is degenerate with the non-magnetic 
{S = 0) state of Bi. This is in agreement with Mielke's 
theoremJl3 according to which the ferromagnetic ground 
state is unique only when a restricted single-particle den- 
sity matrix of the ferromagnetic ground state is irre- 
ducible, and degenerate when it is reducible. That the 
latter case applies for L = 4 can be confirmed by explicit 
calculation of this single-particle density matrix. It is in- 
teresting to note that for A^ = 2 on systems with L > A, 
it can be shown numerically that the ground state is fer- 
romagnetic and unique when the single-particle spectrum 
has two minima. This degeneracy therefore seems to be 
an artifact of the high symmetry of the L — 4 system. 
A^ = 3 : For any ^2 < the ground state is ferromagnetic. 
This is actually the Nagaoka case of one hole in a half- 
filled band. 



For one hole in a half- filled band, the Nagaoka mecha- 
nism leads to a ferromagnetic ground stateJ^J for 1^2 1 — + 0, 
Sigrist and coworkers have-shown that the model is ferro- 
magnetic for all densities J23 and for \t2\ > 0.25, whece the 
band structure has two minima, Miiller-HartmannEZI has 
shown that the low-density limit is ferromagnetic. These 
three limits are indicated in the schematic phase diagram 
shown in Fig. ^. In addition, for 1^2 1 — > 00 the model can 
be mapped onto two decoupled Hubbard chains, which 
cannot be ferromagnetic due to the Lieb-Mattis theorem. 



^ Nagaoka 



Sigrist 



2 Hubbard 
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/ Miiller-Hartmann 



0.25 
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FIG. 4. Schematic U = 00 phase diagram. 



III. VARIATIONAL APPROACH 



The fully polarized state 

i^)= n 4tIo> 



k<ikf. 



is an eigenstatc of Hamiltonian (uh with energy 



^fc 



k<kp 



(4) 



(5) 



where e^ is the single-particle dispersion and kp is the 
Fermi wave vector for the ferromagnetic state. The fer- 
romagnetic state |_F) is certainly unstable if a variational 
state with one flipped spin and a lower energy can be 
found. In order to put good constraints on the extent 
of a ferromagnetic phase, it is important to use as good 
a variational wave function as possible. A particularly 
sophisticated ansatz due to Edwards is defined by 



Af-l 






(6) 



where 



C. Special limits for (7 = 00 



4tW= ^fc.{rn-(.)c'^^^ 



(7) 



For U = 00 and negative ^2, ferromagnetism has an- 
alytically been shown to exist in three different limits. 



creates an up-spin electron in an orbital that is deter- 
mined variationally. The variational parameters are the 



wave vector q and the {N — 1)L amplitudes (pa{i)- For 
i2 = this variational wave function includes the Bethe 
Ansatz as a special case, and is therefore an exact so- 
lution of the simple Hubbard chain.ES This ansatz has 
previously been applied to the ti — ^2 model in order to 
calculate critical U values in Ref. |2^, but here we show 
the results of additional calculations, including the full 
U — oo phase diagram. 

The variaiipnal energy for orthonormal one-particle or- 
bitals readsLJ 



L-l 



E{q,{<f^{£)}) = 

-2ti cos{q) det S^^^ ~ 2^2 cos(2g) det 5"^^^ 



-2ii trS*'^^ - 2i2 trS'f 



2) 






(8) 



(«) 



with overlap matrices S^^I defined as 



{S-S+) = 1 - ^ cos{q6) det S^^'^ 



(5=1 



a, 13 



fa{^)Saf3 ^a{L-5) 



(13) 



For S = 5max = N/2, (y) implies {S- S+) = N, 
whereas for S = S'max — 1, {S^ S^) must vanish. When 
the ferromagnetic state is destabilized, we find that 
(5^5+) « 0.01 — 0.1. That this expectation value is not 
an eigenvalue of S*^ is an indication that the variational 
wave function is not an eigenstate oi H , as opposed to the 
t2 = case. The resulting U = oo phase diagram is dis- 
played in Fig. u. As one can see from the fi gure , all four 
of the analytically treatable limits of Sec. [[IC| are, for 
the most part, reproduced, and the three ferromagnetic 
regions are connected to one another. 



^ij = E<w^/5(^+'^)- 



(9) 



We use the conjugate gradient methodEJ to minimize the 
energy. The derivatives 



dE 



(10) 



can be calculated analytically and are given in Appendix 
^. After each iteration the orbitals are orthonormalized 
using the modified Gram-Schmidt method. The initial 
orbitals, which must be chosen so that they are close 
enough to the final wave functions (pa (^) , are taken to be 
the N — I lowest eigenf unctions of the Hamiltonian dl^) 
with U = and a site impurity at 0. Care has to be taken 
in choosing the particle number N. In particular, the 
corresponding non-interacting ground state must be non- 
degenerate (closed shell) in order to obtain well-behaved 
convergence of the conjugate gradient calculation. 

Here we want to determine the U — oo phase dia- 
gram. Since there is no double occupancy when U = oo, 
faiO) = for all a and the energy gradient with respect 
to (pa{0) can be excluded. We also calculate the total 
spin S = J2i Si of the wave function. Using the commu- 
tator [S'+, S~] — 2Sz and working at a particular Sz, we 
obtain 



N 



N 



1 +(^-5'+). 



(11) 



By applying Wick's theorem to the down spin operators, 
we can write 



{s-s+) = ^J2' 



Aq{t — m) 



i,7n 



a 

(12) 



The elements of the sum can be expressed in term of a 
determinant [cf. Eq. (A12) in Ref. p3; care must be 
taken since there are typographical errors ] , leading to 
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FIG. 5. U — oo phase diagram obtained with ansatz (n) 
from numerical calculations on an L = 100 lattice. The filled 
circles represent stable ferromagnetic states and the open 
squares unstable ones. 

One anomaly is that for small \t2\ the ferromagnetic 
state is destabilized at low density. As we will see, this 
does not occur in the DMRG calculations. This discrep- 
ancy is probably due to the large finite size effects on 
finite lattices at low densities that come from the differ- 
ence between closed and open shells when the system has 
periodic boundary conditions. We have already seen this 
behavior for L = 4 in Sec. HB. The alternating param- 



agnetic and fully polarized states as a function of n near 
the upper boundary of the paramagnetic phase are an- 
other illustration of these finite-size effects. In addition, 
for small n and \t2\ > 8, the ferromagnetic state is un- 
stable in a regime expected to be ferromagnetic by the 
Miiller-Hartmann argument. However, this argument is 
valid in the limit of small n, and here the lowest obtain- 
able n is limited by N/L. 



IV. EXACT DIAGONALIZATION 

We have jPjerformcd nuniericai. diagonalization using 
the LanczosCj and the DavidsonCJ algorithms for chains 
of up to length L = 16, and various numbers of electrons 
A'' and boundary conditions. These methods permit us to 
obtain both the numerically exact energy and the wave 
function of the ground state on a finite cluster. 



A. Boundary conditions 

On small lattices, it is important to analyze carefully 
the effect of boundary conditions. In order to understand 
which boundary conditions should be used, we consider 
the case of TV = 2 electrons on a L = 12 system with 
t2 = —0.1. With periodic boundary conditions the sys- 
tem is never ferr oma gnetic, similar to the four-site model 
discussed in Sec. II B. In fact, this seems to be the case for 



all lattice sizes L. For anti-periodic boundary conditions, 
the single-particle spectrum has 2 degenerate minima at 
Due to Mielke's theoremll3 the model is fully 



± 



12' 



polarized for all C/ > 0. Therefore both boundary con- 
ditions are unable to reproduce the low-density regime 
with a paramagnetic ground state for small U and a fer- 
romagnetic ground state for large U . 

Only with open boundary conditions do we obtain a 
non-magnetic ground state for small U and a ferromag- 
netic one for large U . We find that the critical value 
occurs at Uc ~ 11. Therefore, in order to minimize 
the effect of the boundary conditions, we will take open 
boundary conditions for all exact diagonalization calcu- 
lations as well as for the DMRG calculations described 
subsequently. 



B. Determination of Uc 

In order to determine the critical value of U above 
which the ferromagnetic state is the ground state, we 
start at small U, for which the ground state is not mag- 
netic, and increase U until we reach the fully polarized 
state. The transition point can be determined by exam- 
ining the behavior of the energy Eq. 

As discussed in Sec. [II, the energy of the fully polar- 
ized ferromagnetic state, fi'ferro: does not depend on U 
and is exactly known. Thus, if Eq{U) = Eq{U) = i^ferro 
for all U > U, we identify U with Uc- This can be con- 
firmed by verifying that the lowest eigenvalue is the same 
in all S'z-subspaces, since then the ground state must have 
a degeneracy of 25'max + 1 , or by calculating the expecta- 
tion value of the total spin operator in the ground state, 



(^olS^l^o) = 5I(V'o|S,;S,|V^o) = 5(5- 



1) 



(14) 



For U > Uc one will obtain S — S'max- In Fig. || we 
clearly see that the values for Uc obtained using these 
two criteria are the same. 
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FIG. 6. Eo (open squares) and {ipo\S'^ \tpo) = S'(5'-H) (filled 
circles) for a system of size L = 12 with A'' = 6 and t2 = —0.2. 
The horizontal solid line is the energy of the fully polarized 
state. 



C. Order of the transition 

In order to determine the order of the transition, we 
investigate with very high precision the ground state en- 
ergy Eo{U) around Uc- Since there are many states with 
energy very close to Eq a very large number of iterations 
are needed in the Davidson procedure in order to obtain 
convergence (more than 1000 H\ip) multiplications). 

If the transition is first order, the ground state will 
jump from S = to S = S'max and Eq{U) will have a 
kink at Uc- On the other hand, if the transition is second 
(or higher) order there will be no kink in the energy as 
a function of U and S will smoothly take on all values 
from to S'max- In the thermodynamic limit, a second 
order transition requires that 



,. dEa 

hm — - — 

u^u- dU 



= lim 



dEp 
dU' 



(15) 



i.e. that the derivative of the ground state energy is 
continuous through the transition. 

In order to further clarify this issue we can follow the 
lowest energy state with a particular spin S. Since utiliz- 
ing the S^ quantum number in the exact diagonalization 
program is technically difficult to implement, we follow a 
state of a particular S^ by diagonalizing the augmented 
Hamiltonian 



H' = H + AS^ 



(16) 



in different S^-subspaces with A > 0. For large enough 
A, the lowest energy state within a given Sz sector will 



have the minimum S value. In Fig. ^ we clearly see that 
the spin S of the ground state takes on all intermediate 
values as U is increased. This is an indication that Eq. 
( |l5| ) will be satisfied in the thermodynamic limit and that 
the transition is continuous. Here we have chosen the 
parameters n — 0.5, t2 — —0.2 so that the system is in a 
regime with two Fermi points at U — 0. As will be shown 
in Sec. VI B, the system is a Luttinger liquid for weak U 
at these parameters. The transition is therefore from a 
Luttinger liquid to a ferromagnet. Further evidence that 
the transition is second order based on the behavior of 
the Luttinger-liquid parameters will be given in Sec. M. 
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FIG. 7. Ev,{U) for L = 12, iV = 6 and i2 = -0.2 for 
S = 0, 1, 2, 3. The inset is a blowup of the indicated region. 



D. Small t2 

We have already seen that for ^2 = the fully polarized 
state always has a higher energy than the 5 = state 
unless U — CO. For large but finite [/, we can treat the 
model perturbatively in jj. For ^2 = and n = 1 this 
yields the one-dimensional Heisenberg model. For the 
non-half-filled system and |i2| ^ ii, we obtain, to first 
order in -i. 



H = -ii^(c,}+i^c,:^ + /i 



'cflf 



4i? 
V 



^S.S 



l+l 



7 XI ^^^'cJ+2<T4+l-<TCi+l-<T'Cia' +/i-c.. (17) 



Here A<j<j = +1, \a-a 

t2 



-1 and 



Jpff = rr— 



2tt 



2.2 

— sm nn — sm zttti 



(18) 



is obtained as in Ref. ^. li a = a' the third term of Eq. 
(O) leads to a permutation of the spin part of the wave 
function. The resulting matrix element (proportional to 



effective Heisenberg coupling 



(19) 



where 7 depends only on the filling of the system and, 
in principle, can be calculated. Hence the coupling is 
ferromagnetic when 

-t2>jj (20) 

and for small \t2\ the critical U should behave as 

Uc--\t2\'\ (21) 

A numerical evaluation of Uc in this low \t2\ regime 
(shown in Fig. 0) obtained from exact diagonalization 
is reproduced quite well by this form. 
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FIG. 8. Uc as a function of i2 for L = 10 on a log-log 
scale. The filled circles are for A^ = 4 and the open squares 
for N = 6 while the solid and dashed lines are obtained using 
linear regression. 



V. DENSITY MATRIX RENORMALIZATION 
GROUP 

We now investigate much larger systema2j by apply- 
ing the powerful and already widelyr-^^ed Density-Matrix 
Renormalization Group techniqueJla The DMRG is a 
variational procedure which can be used to obtain the 
energies of the ground state and low-lying excited states 
very accurately, as well as to compute a wide variety of 
equal-time correlation functions. Here we use the finite- 
size algorithm for system sizes of up to L = 140 and keep 
up to 1000 states in the last iteration. The efficiency of 
the algorithm is improved by keeping track of the basis 
transformations in order to calculate a good initial guess, 
for the wave function after adding a site to the system.cfl 
Calculations applied to the ti~t2 model on small systems 



show extremely good agreement with exact diagonahza- 
tion results (up to 10 figures). The total discarded weight 
of the density-matrix eigenvalues provides an estimate of 
the truncation error. In the calculations performed here, 
the discarded density-matrix weight ranges from 10^'* to 
10^^. The estimated error in the DMRG results shown 
here, determined by examining the convergence with the 
number of states kept, is of the order of the plotting sym- 
bol size or less, unless explicitly discussed. 

We have also included the possibility of adding a term 
AS^ to the Hamiltonian. Turning on A > shifts states of 
higher total spin S to higher energies. This shift is known 
for a particular S since S^ commutes with the Hamilto- 
nian. This addition is useful for two reasons. First, since 
the DMRG can only determine a limited number of ex- 
cited states accurately for a given number of states kept, 
it allows more excited states to be accessed within a par- 
ticular S subspace. Second, it allows one to follow states 
of a particular S individually, even if they are not the 
ground state of a particular Sz subspace for A = 0. This 
trick is particularly useful near the ferromagnetic tran- 
sition, where the ground state often has nonzero total 
spin, and states with different total spin are very close in 
energy. The numerical problem of the mixing of different 
total S states when they are near-degenerate in energy is 
also sometimes relieved when A is turned on. 



A. Uc as a function of the density 

In order to determine Uc we investigate the behavior of 
the ground state energy and the expectation value of the 
total spin operator (S^) of the ground state, as in the pre- 
viou sly de scribed exact diagonalization calculations (see 
Sec. [VB ). Here we have chosen to carry out the calcu- 
lations on an L = 50 lattice because the finite-size effects 
are negligible upon further increasing L. In Fig. we 



show results for various values of ^2 as a function of the 
density n. 
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FIG. 9. Uc as a function of density for L = 50 on a 
log- linear scale for three different values of t2 . 



The three curves show three different representative 
behaviors. In all cases, as n — > 1, Uc diverges to reach 
the particular point oi U — oo where all states with dif- 
ferent spin are degenerate. For ^2 = —0.2 (and all other 
cases with one minimum in the single-particle spectrum) 
Uc increases monotonically with n. As n — )■ 0, Uc —^ 0, 
which seems to imply that the problem could be treated 
perturbatively in this limit. However, the relevant pa- 
rameter here is actually Uc divided by n which tends to 
a finite value rather than going to zero at small densities 
(see Fig. PI). 



40 



30 



U,/n 



20 



10 



0.0 



0.2 



0.4 0.6 

n 



0.8 



1.0 



FIG. 10. Renormalized Uc 
L = 50 and i2 = -0.2. 



as a function of density for 



For the two cases t2 = -0.8 and t2 = —2, a local 
minimum appears near n — 0.6 and n = 0.8, respec- 
tively. This minimum in Uc is due to a diverging density 
of states at the Fermi energy of the fully polarized state. 
This occurs at a critical density nc where the Fermi en- 
ergy coincides with the local maximum in e(fc) at /c = 0, 
which exists only for |t2| > 0.25. For ^2 = —2 there ex- 
ists a finite range of densities in which the system never 
becomes ferromagnetic, even aX U — oo. Nevertheless, 
around the critical density there still exists a finite Uc 
with a minimum value at Uc- 



B. Luttinger liquid parameters 

The Luttinger liquid concepta (for a review see Ref. ^ 
or Ref. ^^ is based on a single-particle spectrum with two 
Fermi points {±kp). This is the case in the ii — 12 rnodel 
for all densities when \t2\ < 0.25 and for large enough 
fillings when f2 < 0.25 (see Fig. H). The spectrum is 
linearized in the region of these two points and therefore 
completely specified by the Fermi velocity vp- 

The interaction between electrons can then be writ- 
ten in terms of four scattering processes: backward scat- 
tering (gi), forward scattering ((72), Umklapp scatter- 
ing ((73) connecting the region around +kF with that 



around —kp, and a §4 term connecting states on the 
same branch, either around +fci? or around —kp. The 
Qi term, which is usuaUy neglected since it leads only to 
a renormalization of the Fermi velocity, will be seen to 
be important here. 

We consider a non-half-filled system, where Umklapp 
process can be neglected. Using bosonization, the Hamil- 
tonian can be written in terms of boson field operators 
01/ and their canonically conjugate fields III, {v = p.a for 
charge and spin) . This leads to the spin-charge separated 
Hamiltonian 



H = H„ 



H„ 



2gi 



dxcos 



(22) 



Here a is a short-distance cutoff of the order of the lattice 
spacing and 



H,, ^ dx 



TTUuK^ 2 



2nK^ 



{dx(l>iy) 



(23) 



is the Hamiltonian of an elastic string with eigenmodes 
corresponding to the collective density fiuctuations of the 
fermion liquid. The new, parameters are the charge and 
spin velocities given bycJ 



(24) 
(25) 



and the two coefficients Kp and K„ which determine the 
asymptotic behavior of correlation functions [K^r ~ 1 for 
the Hubbard model due to the SU(2) spin symmetry). 

We can then calculate physical quantities such as the 
spin susceptibility 




2K, 

TT U„ 



or the density-density correlation function 



(26) 



{5n{x)5n[0)) 



Kn 



(ttx)^ 



Ai cos{2kFx)- 



;l+^p log3/2 : 



(27) 



where dn{x) — n{x) — {n{x)). 



1. Determination of Kp 

To obtain Kp using DMRG we compute the Fourier 
transform 



^('?) = ;^E^"''^-c(^) 



(28) 



of the charge-charge correlation function ( p7| ) and then 
take the first derivative at q = 



dC{q) 



K, 



(29) 



?=o 



The q = {) derivative is proportional to the coefficient of 
the X/x"^ term in Eq. ([27|). This method of extracting 
Kp has been shown to yield accurate results for the ex- 
actly solvable case ^2 = and should be valijijas long as 
the system is in the Luttinger liquid phase O Since we 
work with open boundary conditions we have to use zero 
padding for \i\ > L and average the correlation function 
to reduce boundary effectalj: 



A^ave(^) = — V N{io + m,io+i + m) 

n _ -^ — ^ 



rn— 



where ig -I- i/2 = L/2 and 



NiiJ) = (riinj) - (ni)(nj). 



(30) 



(31) 



The quantity ria is taken to be large enough so that 
-^ave(^) does not depend strongly on io and Ua', typically 
we take ria « 6. 

The correlation functions for a system of size L — 100 
and ^2 = —0.2 with density n = 0.5 are given in Fig. |lj. 
One can see that the slope at g = is well-behaved and 
decreases monotonically with increasing U. 
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' FIG. 11. Fourier transform C(g) of the density-density cor- 
relation function for systems of size L = 100, t2 ~ —0.2, 
n = 0.5 and t/ = 0, 1,...,8. 

Therefore Kp can be accurately calculated using Eq. 
(P9|). In Fig. |2| the numerical results for ^2 = —0.2 
are compared with the well-kspwn values for the sim- 
ple (^2 = 0) Hubbard chain,L3 where the appropriate 
parameter in the weak-coupling regime is U/vp. It is 
seen that the two curves agree very well for weak cou- 
plings, although the Fermi velocities differ appreciably: 
vp{t2 = 0)/vp{t2 = -0.2) « 2.3. The deviation be- 
tween the two curves increases for larger couplings. For 
t2 = 0, Kp goes asymptotically to the value i as C/ — > 00, 
whereas for ^2 — —0.2 Kp reaches i at a finite U whose 
value agrees quite well with the Uc calculated above. 
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FIG. 12. Kp as function of t/ forn = 0.5. The filled circles 
are for t2 — —0.2 and the dashed line is for t2 — 0. 



In Fig. ^, we show the finite-size scaling of the spin 
gap for a system with density n = 0.5 and ^2 = —0.2 for 
increasing values of U (from 1 to 5). For these parame- 
ters, Uc « 7.55. The excitations are gapless, as expected, 
and the spin velocity u^ decreases to with increasing 
U. Here K^r — 1 due to the SU(2) invariance, and the 
susceptibility ( |2q ) is then directly proportional to u~^. 
Hence, as can be seen in Fig. |lj, x diverges when ap- 
proaching Uc- A diverging susceptibility is an indication 
of a second-order transition, in accordance with the anal- 



ysis of Sec. IV C . The critical exponent 7 is defined near 



the transition by 



X' 



(35) 



where u = \U — Uc\- We obtain 7 = 2.0 ± 0.1 by fitting 
the results shown in Fig. y^, where the error is from the 
least-squares fit. 



2. Spin and charge velocities 

To determine the velocity of the spin and charge ex- 
citations, we use finite-size scaling of the corresponding 
energy gap 

TT 



L+1 



(32) 



where Afc is the finite interval between two adjacent k- 
points forj-open boundary conditions. The spin gap is 
defined aa£d 



A, = EoiS ^ 1) - EoiS ^ 0) 



(33) 



where Eq(S) is the lowest eigenvalue for a system with 
spin S, and the charge gap is given by 

Ap = i [EoiN + 2) + E„iN-2)- 2Eo{N)] (34) 

where Eo{N) is the ground state energy for N particles 
and Sz = 0. 
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FIG. 13. Spin gap for different values of U for systems with 
fa = -0.2 and n = 0.5. 
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FIG. 14. Spin susceptibility as a function of U for a system 
with t2 = —0.2 and n = 0.5. 

We have also calculated the charge velocity. In Fig. 
p^ we show the results for the same parameters as above 
and compare them with the case ^2 = (obtained fol- 
lowing Ref. KO). We notice that for ^2 — —0.2 the charge 
velocity, Up, is strongly renormalized from the t2 = 
value as U is increased, even when U is rescaled by the 
Fermi velocity. Eq. (|2^) suggests that this behavior could 
be governed by the g^ term which describes scattering 
between states on the same branch of the spectrum and 
whose effect is usually taken to be unimportant. The 
increased importance of the 174 interaction in a system 
with a tendency towards ferromagnetism is understand- 
able since the system's response to an external magnetic 
field is described in terms of the operator 



= 



E-[i. 



-k)cr^{kf 



hk}a 



-{-kF + 



k)a% 



kF-\-k)a 



k,a 



(36) 



which involves states on the same branch of the spectrum. 
This indicates that (74 must appear in the renormalization 



10 



of the corresponding response function. Hence the 54 
term must be relevant for the ferromagnetic fixed point, 
and we would therefore expect u^ to be strongly rescaled. 



of the order of the system size as the transition is ap- 
proached, so that we can at most say that the best fit to 
the data occurs with a relatively small value of z/, v k, 0.2. 
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FIG. 15. The charge velocity Up as a function of U for 
^2 = (full line) and ^2 = -0.2 (filled circles). 



C. Spin-spin correlation function 

We can also study the behavior of the spin-spin corre- 
lation function near the transition. In order to minimize 
effects of the open boundary conditions, we average over 
a number of sites for a given distance, as was done for 
the density-density correlation function in Sec. VBl, 
and consider 



'S'ave(^) = V" S{io 



m,io 



m=0 



where io + i/2 = L/2, 



4 and 



Siz,j)^{SrS+). 



(37) 



(38) 



Fig. nfl shows the result for values of U near Uc- We find 
that S'ave(^) is positive definite, indicating ferromagnetic 
correlations, and can be well-fitted by the form e^^^^, 
with ^ diverging as the transition is approached. This 
is seen by the linear behavior of S'ave(^) on the semi- log 
plot, with decreasing slope as U increases and the tran- 
sition is approached. We can define the critical exponent 
associated with the divergence of the correlation length 
as 



c 



(39) 



where u = \U — Uc\. Unfortunately, near the transition 
mixing of energetically close states make it numerically 
difficult to accurately calculate Sa,ve{^), even using 7 it- 
erations and keeping up to 800 states. There tend to 
be systematic errors, which we minimize by limiting the 
system size to i = 40. However, ^ then quickly becomes 



FIG. 16. The averaged spin-spin correlation function 
Save(^) for a system with L = 40, t2 ~ —0.2 and n — 0.5 for 
U — 7.60, 7.65, . . . , 8.00 and 10, plotted on a semilog scale. 

The total spin S of the ground state is also related to 
this correlation function (sec Eq. dlj)) via 

(S^) = SiS + 1) = 5.(5. + 1) + ^(5r5+), (40) 

which for a translationally invariant system becomes 

S{S+l)^S,{S, + l) + LY,S{i), (41) 

where S{i) = {S^ S^^) is independent of i. Therefore, 
the ferromagnetic order parameter is 



S/L ^ L-1/2 



E^w 



1/2 



(42) 



If J2i S{£) is finite for L ^ 00, the system is disordered; if 
it is proportional to L, there is long-range ferromagnetic 
order; and if it follows a power law L'^~'^ with 1 < rj < 2, 
the system is at the critical point with critical exponent 
77 and S{£) ~ fi~^ . In the ordered ferromagnetic phase, 
we do find a nonzero value of S{£) at large distances as 
seen in Fig. |l^, consistent with this picture, but we have 
not been able to determine the critical exponent rj from 
numerical calculations because of the poor convergence 
of the DMRG at the critical point and because of uncer- 
tainty in Uc- 



VI. PHASE DIAGRAM 
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A. [/ = CO phase diagram 

The phase diagram for U — oo can be determined us- 
ing DMRG calculations. While one could, in principle, 
start with a Hilbert space in which double occupancy is 
explicitly excluded, here we simply set U = 10^ to mimic 
the infinite-C/ limit and find that the accuracy is quite 
good since the DMRG integrates out high energy scales 
automatically. In order to decide if a point in the n — t2 
plane is ferromagnetic or not, we compare its DMRG 
energy Ed with the known ferromagnetic energy -Bferro 
and, in addition, calculate (S^) = S{S + 1) and com- 
pare it with 5'max- If 5 « 5'max and Ed > £^ferro for 
given ^2 and n we conclude that the system is ferromag- 
netic for these parameters, whereas if we find an energy 
Ed lower than -Eferro and 5 « 0, we conclude that it 
is non-magnetic. Near the boundary between these two 
regions we sometimes find partly polarized states, i.e. 
Ed < £^ferro but S > 0. It is difficult to determine 
the nature of these states because there are two possible 
causes for the partially polarized value of S. One pos- 
sibility is that the system undergoes a continuous phase 
transition as a function of n or t2 at U ^ oo. However, 
near a ferromagnetic phase transition, near-degeneracy of 
states leads to a mixing of states in the diagonalization 
step of the DMRG procedure. Therefore, the presence of 
partially polarized states could also be due to numerical 
effects. This mixing can also be seen in that the values 
of S obtained in general take on continuous values that 
lie between the discrete values of S allowed by the finite 
number of electrons in the system. These problems are 
similar to those that occur near th e ferr omagnetic tran- 
sition at finite U discussed in Sec. IVC, where we used 



extremely accurate exact diagonalization calculations in- 
stead of DMRG to solve them. 



romagnetism. The three limiting cases of Sec. 11 C are 
reproduced and the corresponding ferromagnetic regions 
are, in fact, connected. We also notice that for sufficiently 
large \t2\, the system is not magnetic at intermediate den- 
sities. This extends the limit ^2 —^ — oo in which the 
system consists of two uncoupled Hubbard chains which, 
must be non-magnetic due to the Lieb-Mattis theorem,t2l 
to a finite region. In this region, the system behaves effec- 
tively like an uncoupled two-chain model. We will later 
present evidence that this region in which the system be- 
haves as two uncoupled Hubbard chains extends to the 
low-U phase diagram. Note that this phase diagram is 
qualitatively similar to that obtained using the Edwards 
ansatz (see Fig. ||), except that here the entire ^2 axis is 
ferromagnetic for small densities. 



B. Low- (7 phase diagram 

We now turn to the general question of the low U 
phase diagpam within a weak-coupling analysis. Balents 
and FisherES have analyzed the weak-coupling phase di- 
agram of the two-chain Hubbard model using RG and 
bosonization. Their calculation is generic for a system 
with four Fermi points. They obtain coupled RG equa- 
tions and integrate them numerically to find the different 
fixed points, which they then analyze using bosonization. 
The possible phases can be classified by the number of 
charge and spin modes which are gapless. A phase with 
a gapless charge modes and (3 gapless spin modes is de- 
noted CaS/3, where a and /3 can take on integer values 
from to 2. An adaption of their two-chain phase dia- 
gram to the ti-t2 model in the t2-n plane is shown in Fig. 
nH, with the ^2 < 0, n < 1 quadrant shown in more detail 
in Fig. H. 
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FIG. 17. Phase diagram for U — oo and L = 30 ob- 
tained with DMRG. The filled circles indicate a fully polarized 
ground state, the open square a non-magnetic ground state 
and the stars a partially polarized ground state. 

As seen in Fig. |l^, we obtain a large region of fer- 




FIG. 18. Low-U phase diagram obtained by adapting the 
results of Ref. 32. Some small regions of additional phases 
near the solid lines are not shown. 



As discussed in Sec. VB, a Luttinger liquid phase 
(ClSl) is expected for small \t2\ in the region where 
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the U = system has two Fermi points. At half-fiUing 
2kF = TT and umklapp processes cause the system to 
open a charge gap, so that the phase is that of the one- 
dimensional Heisenberg model (COSl). 
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FIG. 19. Low-U phase diagram restricted to negative t2 
and n < 1. The filled circles are for parameters for which no 
spin gap has been found with DMRG at f/ = 2, while the 
open circles are for parameters where a spin gap has been 
found. 

When the Fermi surface has four points, namely ±fci?i 
and ±fci?2, the effective low-energy model has four lin- 
earized single-particle branches ei (k) ^ T'>^Fi{k — kp-^) 
and e2(fc) = iw^aCfc — kp^) (see Fig. |20| ). 
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FIG. 20. Efi'ective low-energy model for a Fermi surface 
with four points. 

This is equivalent to the two-band model of Ref. ^2 for 
which there are also four Fermi points, namely ifcfT for 
the bonding, and ±kp for the antibonding band. The 
only difference here is that the inner bands, which are 
denoted ei{k) in the present case and originate from the 
antibonding band in the two-chain model, have opposite 
velocities. The correct mapping of the Fermi points be- 
tween the two models is therefore 



±kp2 



Tk% 
±kl. 



Provided that this mapping is performed, the perturba- 
tion expansion of the ti — ^2 model and the two-chain 



Hubbard model is exactly the aapie at low energy, as al- 
ready pointed out by Fabrizio.E2l Therefore, we can sim- 
t adapt the results from Ref. |2| to our case. In Fig. 
and Fig. |lj, the thick solid line represents the critical 
density ric for which the Fermi surface splits into four 
points (cf. Fig. ||). Exactly on this line, Balcnts and 
Fisher predict a CISO phase. For slightly smaller densi- 
ties, vp-^ is much smaller than vp^ leading first to a C2S2, 
then to a C1 S2 p hase. (These three phases are not de- 
picted in Fig. |l^ since they have a small extent.) When 
Vp^ is comparable to vp^ the weak-coupling RG leads to a 
large region of a CISO phase. This phase is a doped spin- 
liquid phase with a spin gap in-which power-law pairing 
and CDW correlations coexist .E3 When 2k p^ — n within 
this region, indicated by a dashed line, umklapp pro- 
cesses in the corresponding bonding band can open a 
charge gap. For the two-chain model, Balents and Fisher 
predict a C1S2 phase along some of this line, with the 
result being sensitive to the initial conditions in the RG 
equations, i.e. on vp and the initial couplings. 

In order to investigate the validity of this weak- 
coupling phase diagram, we have calculated the spin and 
charge gaps using the DMRG for different system sizes 
at small U (we choose U = 2) and a number of t2 and n 
values. Due to the weak coupling and small size of the 
gaps, very high precision is necessary in the DMRG pro- 
cedure. We use up to 8 finite-size iterations and keep up 
to 800 states in the last iteration. The presence of a spin 
gap in the extrapolated L — > oo limit is indicated by an 
open circle in Fig. 19 and the absence of a spin gap by a 



solid circle. 

Figs. |2l| and g2| show the finite-size scaling of the charge 
and spin gap. The filled circles represent the spin gap 
and the open squares the charge gap while the lines are 
quadratic regression in (i) between the points. In gen- 
eral, the finite-size corrections to the spin and charge gaps 
for a system with open boundaries can be represented as 
a power series in 1/L. When a gap is present, the domi- 
nant correction is usually l/L^,|-and when the system is 
gapless, the dominant correction^ is lAL. While this be- 
havior is generally seen in Figs. ^ and |2^, there are cases 
with a small gap with an obvious positive quadratic cor- 
rection, but also a substantial linear term. In addition, 
there is scatter of up to the order of the symbol size in 
some of the curves which we believe is due to additional 
finite-size effects which can oscillate with the number of 
particles. 
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FIG. 21. Finite-size scaling of charge (open squares) and spin (filled circles) gap for different parameters. 



For ^2 = —0. 2,71 = 0.5, the clearly vanishing charge and spin gaps confirm that the system behaves as a Luttinger 
Uquid (see Sec. VB). For parameters close to the line of critical n, namely t2 — —0.5, n = 0.4 and i2 = —0.5, n = 0.8, 
we find no spin and charge gap, in contradiction with the adapted weak-coupling phase diagram. However, a finite 
U could renormalize the band structure, leading to a shifting of the line o£-transition from two Fermi points to 
four. Such a shift is seen with increasing U in the two-chain Hubbard modelEj Precisely on this transition line, for 
t2 — —0.38, n = 0.8, a spin gap would be predicted, but we do not find one. However, this phase might be hard to see 
numerically since it occurs only exactly on the line, or nught not be present for finite U. Such a phase has also not 
been found numerically in the two-chain Hubbard modcl.CJ 

For the case t2 = —1.25, n = 0.8 we find a spin gap, as predicted. For ^2 ~ —0.5, n = 0.2, where a spin gap is also 
predicted, we have taken U = 0.5 rather than U ~ 2 because the ground state is not paramagnetic at the larger 
U value due to the proximity of the ferromagnetic transition. At the smaller U, finite size effects and numerical 
problems make it difficult to definitely determine whether or not there is a spin gap. Nevertheless, since there is a 
strong quadratic correction, we conclude that a spin gap is probably present. 

Near the line 2kp2 = tt, (i.e. at ^2 — — 1,« — 0.4) the system has a small spin gap, consistent with the CISO phase 
but not a possible C1S2 phase. However, as pointed out above, the existence of the C1S2 phase in the weak-coupling 
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calculation is dependent on the initial conditions, so it .aaay not be present along all of the dashed line. In the 
two-chain Hubbard model, numerical DMRG calculationsc3 do find evidence for this phase for some fillings, even at 
intermediate to strong U . More work would have to be done to determine in detail some of the finer structure of the 
phase diagram, but the difficulty of the calculations and the finite-size scaling preclude a more detailed investigation 
here. 

For ^2 = —2.0 and n — 0.4, we find that both the spin and the charge gaps vanish, in contradiction to the weak- 
coupling phase diagram, which would predict a CISO phase. Thi s occu rs in the regime which is paramagnetic at large 
U because the system behaves as two uncoupled chains (see Sec. VIA ). We therefore suspect that the strong-coupling 
behavior extends to weak coupling, and that the system is in a 2 x ClSl = C2S2 phase here. The phase boundary 
of the paramagnetic strong-coupling phase, indicated by the solid diamonds, is also sketched in. It remains to be 
determined how much this phase boundary changes in going from strong to weak U . 
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FIG. 22. Finite-size scaling of charge (open squares) and spin (filled circles) gap for difi'erent parameters. 
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It should also be noted that Kuroki et al.Ell have stud- 
ied the ti-t2 model numerically using Projector Quantum 
Monte Carlo and the DMRG at half-filling and ta = -0.8 
and find a spin-gapped metallic phase with dominant 
pairing correlations at weak U, in agreement with Fig. 
|l8| and Ref. ^. As U is increased they find a transi- 
tion to a spin-gapped insulating phase also in agreement 
with Ref. ^. Here we have concentrated on n < 1 in 
the numerical work since the ferromagnetic phase is not 
present at n = 1. In very recent work that we became 
aware-of as this manuscript was being completed, Arita 
et al.Ea, have calculated the spin gap at n = 1.0 and 
n = 0.5 and ij = -0.55, -0.8, -2.0 using the DMRG. Of 
particular interest for the weak-coupling phase diagram 
are the n = 0.5, ^2 = —0.8, —2.0 points which arc in a 
paramagnetic phase. For the ^2 = —0.8 point, they find 
a very small but finite gap at f7 = 8, consistent with 
a CISO phase. However, we believe that their data are 
also consistent with a vanishing spin gap, which would be 
consistent with the C1S2 phase present along the dashed 
line in Fig. |l^. At ^2 = —2.0 and U = 16, they find a van- 
ishing spin gap, in contradiction with the weak-coupling 
prediction of a gapped CISO phase at these parameters, 
but consistent with the 2 x ClSl phase proposed above. 



VII. CONCLUSION 

We have studied an extended version of the conven- 
tional one-dimensional Hubbard model in order to inves- 
tigate the mechanism for ferromagnetism in an itinerant 
electron model. The added term, which involves hopping 
between next-nearest neighbor sitesp— precludes the ap- 
plication of the Lieb-Mattis theoremE3 which excludes a 
ferromagnetic ground state in the one-dimensional Hub- 
bard model. Indeed, we do find a ferromagnetic phase 
in a wide region of parameters at large enough U in the 
regime with ia negative and n < 1 (which is equivalent 
to the ^2 > 0, n > 1 region through a particle-hole trans- 
formation) . 

Using exact diagonalization, the Density-Matrix 
Rcnormalization Group and the Edwards variational 
ansatz, we have shown that the three different mecha- 
nisms for ferromagnetism obtained by taking special lim- 
its at t/ = oo (the Nagaoka state, the limit of vanishing 
density and the limit of very small ta) are all connected 
in the same phase in the t2-n plane. For large negative 
ta/ii , there is a paramagnetic region in the large-C/ phase 
diagram in which the system behaves like two indepen- 
dent Hubbard chains. This region extends to i2 ~ —1.5 
at intermediate n. 

The critical interaction strength for the ferromagnetic 
transition, Uc, has three characteristic behaviors as a 
function of n. When > ^2 > —0.25, there is one min- 
imum in the single-particle dispersion, e(fc), and Uc in- 
creases monotonically with n. For —0.25 > jial > —1.5, 
there are two minima in e(fc), and there is a local min- 



imum in Uc when the Fermi level of the fully polarized 
ferromagnetic state is at the singularity in the density 
of states corresponding to the local maximum in e(fc). 
Finally, when \t2\ < —1.5, there is an intermediate re- 
gion on n for which there is no ferromagnetism, even at 
U = oo, but there is a finite local minimum in Uc when 
the Fermi energy is at the local maximum in e(fc). 

This leads us to the question of what general properties 
are required in order to obtain metallic ferromagnetism 
in this model. The general picture is that with hole dop- 
ing, ^2 rnust be less than zero in order to obtain a fer- 
romagnetic state. When this condition is satisfied, the 
ferromagnetic state occurs over a wide range of parame- 
ters, with, in some cases, quite small Uc- The mechanism 
for ferromagnetism can be motivated from a local point 
of view, in that when ta is negative the triangular struc- 
ture of the chain frustrates the antiferromagnetic order 
(a generic effect for lattice models) . That this frustration 
can lead to a ferromagnetic ground state can be seen on 
small cluster calculations. 

Another point of view emphasizes the importance of 
thoJprm of the single-particle density of states. Wahle et 
al.J3 for example, emphasize that a necessary condition 
for ferromagnetism is an asymmetric density of states, 
with a strong singularity and a larger density of states 
in the lower part of the band. Our results here also sup- 
port these ideas. For ^2 = 0, the density of states is 
symmetric and there is no ferromagnetism. When ^2 < 
the density of states becomes asymmetric, with the e^^/^ 
Van Hove singularity at the lower band edge gaining in 
weight. For ^2 < —0.25, the presence of the additional 
Van Hove singularity at the ferromagnetic Fermi level 
further stabilizes the ferromagnetic ground state. 

The weak-coupling behavior of this model has also 
proven to be quite interesting. For weak negative t2, the 
low-energy effective behavior of the model does not differ 
qualitatively from that of the one-dimensional Hubbard 
model. For weak U, we do indeed find that the model is 
well-described as a Luttinger liquid, and have been able 
to extract the Luttinger- liquid parameters, Kp, Up and Mo- 
using the DMRG. However, unlike the one-dimensional 
Hubbard model, the Luttinger liquid state of the ti-t2 
model undergoes a transition to a ferromagnetic state at 
finite U . While it is clear that the breakdown of the 
Luttinger liquid is not described within the usual weak- 
coupling picture, we have tried to indicate how the break- 
down occurs within this picture. The spin velocity, which 
goes to zero asymptotically as C/ — > oo in the t2 — Q case, 
becomes zero at finite Uc for ^2 < 0. The ferromagnetic 
susceptibility, which is inversely proportional to the spin 
velocity in a Luttinger liquid thus diverges at the tran- 
sition, implying that the transition is second order. By 



Uc 



we 



fitting this susceptibility with the form \U 
obtain a critical exponent 7 = 2.0 ± 0.1. 

In addition, we have calculated the spin-spin correla- 
tion function in the vicinity of the transition and find that 
it becomes ferromagnetic and exponentially decaying just 
below the transition. The correlation length grows as the 
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transition is approached from below, which is consistent 
with a second order transition. We have attempted to 
extract a critical exponent v by fitting the correlation 
length to a form ^ ^ \U — Uc\~^ , but find that it is diffi- 
cult to extract an exponent due to convergence problems 
which limit the maximum lattice size near the transition. 

An examination of the behavior of the ground state 
energy as function of U near the transition using exact 
diagonalization suggests that Eq{U) becomes smooth in 
the thermodynamic limit and provides further evidence 
that the transition is second order. 

Finally, we have investigated the very rich low-C/ phase 
diagram. For a large region of parameters the low-t/ 
phase diagram of the ii — ^2 model can be mapped to that 
of the two-chain Hubbard model. For sufficiently large 
I f 2 1 and a wide range of n < 1 , we confirm numerically the 
existence of the doped spin-liquid phase (CISO) predicted 
by weak coupling RG, which is the one-dimensional ana- 
log of a superconducting phase. In addition, we have 
presented evidence for the existence of a new 2 x ClSl 
phase (not found in the weak coupling treatment) in a 
region in which we think that the ii — ^2 model behaves 
as two uncoupled Hubbard chains. Because we have also 
found some additional discrepancies between the weak- 
coupling phase diagram and the numerical calculations, 
more work needs to be done to clarify the details of the 
phase diagram for the lattice model. 



Ai: 


-2t2 + v/4t2 + 8*2 


E: -2ti 




-2t2 - ^4.tl + Stf 


-2h 
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4t2 
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-4t2 


2ti 







2h 



TABLE II. Eigenvalues of the square model for N = 2. 
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E : -2t2 
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-^ztl + ul 
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TABLE III. Eigenvalues of the square model for A'^ = 3. 
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-1 


Bi 


1 


-1 


1 


1 


-1 


B2 


1 


-1 


1 


-1 


1 


E 


2 





-2 









TABLE I. Character table of the group C4v 
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APPENDIX A: EDWARDS ANSATZ 



The gradient of E with respect to the one-particle orbitals faij) ~ Vaj can for real wave functions be simplified to 

= 2{x\{H-E)^-—\x). (Al) 



This expression can be evaluated using Wick's theorem, yielding 

dE 



dLpaj "J "^ "^ 



= F'^,+FL+K (A2) 



where 



^4- = -tiCos(<7)det5« J ^^;3j+l[5W"]^„ + ^^,_i[5(i)"']„/3 - 2(^„, I 

-t2COs(2<j)det5(2) J ^^^J.+2[S(')"']/3a+^/3J-2[^''^"\/3-2(^„J 
-*2 C0s(2g) < ^ (^/3j [S^al + Sal) ~ "^"^+2 " Vaj-2 



Faj = C^'/'aO \ 5jO - ^ <^/30'y3/3j 



with overlap matrices S^L defined as 



sl:l = J2v*MMe + ^)- (A3) 



e=i 
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